Demonstration

Jean-Romain - in partnership with the Quebec ministry - has been developing approaches to detect, measure, and evaluate road detection from lidar data. His approach leverages the lidR package, which has seen massive uptake and use globally. lidR, also developed by Jean-Romain, is free and open source and is considered a key resource for lidar data processing globally.

Algorithm overview

The implementation of this algorithm is a work in progress. Function names are likely to change, methods are likely to evolve, and debug options are likely to be removed. Simply described, the current road detection and measurement algorithm measure_road() leverages the lidR processing engine to extract relevant segments of lidar data from a point-cloud with a paired linear road network. The approach is in ongoing development, but follows the following steps:

  1. A lidar dataset and road segment and provided

  2. The road segment is buffed and used to clip relevant lidar data

  3. The road is divided into incremental segments

  4. Structural lidar metrics are calculated at each segment used to define the location with the highest probability for the road

  5. Important road measurements are calculated including drivable width and embankment width.

  6. If the algorithm is confident that the road exists, a new (corrected) road segment is produced.

Example road

This example uses a section of SPL data and road network shapefile from the Nipissing Forest.

#--- load in required libraries ---#
library(lidR)
library(sf)
library(MFFProads)
library(dplyr)
library(leaflet)

#--- read in laz and road shapefile ---#
infolder <- "I:/Example_data/LAZ_files"
ctg <- readLAScatalog(infolder, filter = "-drop_withheld -keep_random_fraction 0.25")
roads <- st_read("E:/SPL_Roads/data/Nipissing/tst_roads_clip3.shp")
roads = st_cast(roads, "LINESTRING")
mapview::mapview(list(ctg,roads), layer.name = c("Catalog","Roads"),
                 color = "red", map.type = "Esri.WorldImagery")

Detect road and calculate metrics

Individual road segments can be visualized alongside co-located SPL data.

#--- test on single road ---#
road <- roads %>%
  dplyr::filter(FID %in% c(10))

Detecting the road

We define that we want to find the road using options(). An example of how data are visualized is provided below. Note the different metrics used. A spike in the purple lines denotes an increased probability of road for a given point in the transect.

#--- find the roads ---#
param = mffproads_default_parameters
param$extraction$road_buff = 25

corrected_road <- MFFProads::measure_road(road = road, ctg = ctg, dtm = dtm, water = water, param = param)
mapview::mapview(list(road, corrected_road),
                 layer.name = c("Inaccurate", "Corrected"),
                 color = c("red", "blue"), map.type = "Esri.WorldImagery")

Example of road detection

Measuring the road

Following detection of the road the algorithm can also measure important road parameters.

Example of road measurement

Derive road metrics

Metrics about the road can also be measured, giving details about drivable quality and embankments. A final STATE of the roads is provided with the following gradient:

  1. Road exists, very confident
  2. Road exists, some confusion
  3. Road doesnt exist, some confusion
  4. Road doesnt exist

Example of scoring

Example using multiple roads

This example uses a section of SPL data and road network shapefile from the Nipissing Forest.

Road segments can be visualized alongside co-located SPL data.

roads <- st_read("E:/SPL_Roads/data/Nipissing/tst_roads_clip3.shp")
roads = st_cast(roads, "LINESTRING")
#--- test on single road ---#
roads <- roads %>%
  dplyr::filter(FID %in% c(1:120))

Detecting the roads

We define that we want to find the roads using options(). An example of how data are visualized is provided below.

#--- find the roads ---#
param = mffproads_default_parameters
param$extraction$road_buff = 25

corrected_roads <- MFFProads::measure_roads(road = roads, ctg = ctg, dtm = dtm, water = water, param = param)
mapview::mapview(list(roads, corrected_roads),
                 layer.name = c("Inaccurate", "Corrected"),
                 color = c("red", "blue"), map.type = "Esri.WorldImagery")